#----- Load posterior samples
load("result_Effects")

#----- Number of posterior samples
len <- dim(result)[1]

#----- Min, 2.5%, Mean, 97.5%, Max of posterior samples
x1 <- apply(result[,1:9], 2, min)
x2 <- apply(result[,1:9], 2, function(x) sort(x)[len*0.025])
x3 <- apply(result[,1:9], 2, mean)
x4 <- apply(result[,1:9], 2, function(x) sort(x)[len*0.975])
x5 <- apply(result[,1:9], 2, max)

TE <- c(x1[1],x2[1],x3[1],x4[1],x5[1])
JNIE123 <- c(x1[2],x2[2],x3[2],x4[2],x5[2])
NDE <- c(x1[3],x2[3],x3[3],x4[3],x5[3])
NIE1 <- c(x1[4],x2[4],x3[4],x4[4],x5[4])
NIE2 <- c(x1[5],x2[5],x3[5],x4[5],x5[5])
NIE3 <- c(x1[6],x2[6],x3[6],x4[6],x5[6])
JNIE12 <- c(x1[7],x2[7],x3[7],x4[7],x5[7])
JNIE23 <- c(x1[8],x2[8],x3[8],x4[8],x5[8])
JNIE13 <- c(x1[9],x2[9],x3[9],x4[9],x5[9])

#----- Box plots of the effects
pdf("BoxPlot.pdf", height=4, width=9)
boxplot(list(TE,JNIE123,NDE,NIE1,NIE2,NIE3,JNIE12,JNIE23,JNIE13),range=0, las = 1, names = expression(TE[],JNIE[123],NDE[],NIE[1],NIE[2],NIE[3],JNIE[12],JNIE[23],JNIE[13]))
stripchart(list(result[,1],result[,2],result[,3],result[,4],result[,5],result[,6],result[,7],result[,8],result[,9]), vertical = TRUE, method = "jitter", pch = 21, col = c(rgb(1,0,0,alpha=0.1),rgb(1,0,0,alpha=0.1),rgb(1,0,0,alpha=0.1),rgb(0,0,1,alpha=0.1),rgb(0,0,1,alpha=0.1),rgb(0,0,1,alpha=0.1),rgb(0,1,0,alpha=0.1),rgb(0,1,0,alpha=0.1),rgb(0,1,0,alpha=0.1)), add = TRUE)
dev.off()

